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ABSTRACT 

This work introduces a new fitting formalism for isophotes which enables more accurate modelling 
of galaxies with non-elliptical shapes, such as disk galaxies viewed edge-on or galaxies with X- 
shaped/peanut bulges. Within this scheme, the angular parameter which defines quasi-elliptical 
isophotes is transformed from the commonly used, but inappropriate, polar co-ordinate to the 
‘eccentric anomaly’. This provides a superior description of deviations from ellipticity, better cap¬ 
turing the true isophotal shape. Eurthermore, this makes it possible to accurately recover both the 
surface brightness profile, using the correct azimuthally-averaged isophote, and the two-dimensional 
model of any galaxy: the hitherto ubiquitous, but artificial, cross-like features in residual images 
are completely removed. The formalism has been implemented into the IRAE tasks Ellipse and 
Bmodel to create the new tasks ‘Isofit’, and ‘Cmodel’. The new tools are demonstrated here 
with application to five galaxies, chosen to be representative case-studies for several areas where 
this technique makes it possible to gain new scientific insight. Specifically: properly quantifying 
boxy/disky isophotes via the fourth harmonic order in edge-on galaxies, quantifying X-shaped/peanut 
bulges, higher-order Eourier moments for modelling bars in disks, and complex isophote shapes. 
Higher order (n > 4) harmonics now become meaningful and may correlate with structural 
properties, as boxyness/diskyness is known to do. This work also illustrates how the accurate 
construction, and subtraction, of a model from a galaxy image facilitates the identification and 
recovery of over-lapping sources such as globular clusters and the optical counterparts of X-ray sources. 

Keywords: galaxies: individual (ESO 243-49, NGC 2549, NGG 3610, NGG 936, NGG 4111) - galaxies: 

fundamental parameters - galaxies: photometry - galaxies: structure - techniques: image 
processing 


L INTRODUGTION 

Galaxies come in a wide variety of shapes and sizes. 
Among the numerous efforts throughout the years to 
impose some order among their ever growing numbers, 
one particular classification scheme remains persistent 
to this day: grouping galaxies by their shape as it is 
seen at optical and infrared wavelengths. This is es¬ 
sentially a structural classification - while some galax- 


ies, known as 

‘late ty 

pe’ on the Hubble/Jeans sequence 

(Jeans 1919 

IHubble 

1926, Jeans 

1928), display disks. 

spiral arms. 

oars anc 

bulges, others are more relaxed 


systems, and are commonly referred to as ‘early type’, 
elliptical or lenticular galaxies. However, even ‘elliptical’ 
galaxies are rarely truly elliptical. Their isophote shapes 
often deviate from pure ellipses, in a characteristic way. 
These deviations, the most common being referred to 
as ‘boxyness’ or ‘diskyness’, originate from the structure 
of the stellar orbits that make up the galaxy. Because 
of the physical link between isophote shape and galaxy 
(structural) properties, quantifying these deviations pro¬ 
vides a valuable tool to study galaxies. Numerous works 
have revealed correlations between boxyn ess/diskyness 


and physical properties such as kinematics (Garter 

1978 

Davies et al.|1983|Lauer|I9851|Carter|I987,|J3ender 

rm 

Peletier et al.|| 

11990, Jaft'e et al.||I994p, brig] 

titness prohles 
d X-ray prop- 

(IlNieto et al.|1991D and even global radio an 


erties ( [Bender et al.||I989 ). 

However, the fornialisni used in the past (and to this 
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day) to describe boxyness/diskyness in isophotes has lim¬ 
ited applicability. In particular, it is only efficient when 
applied to galaxy isophotes which are relatively well de¬ 
scribed by pure ellipses and the deviations from elliptic¬ 
ity are small (a few per cent). This has limited isophotal 
structure studies to only a subsample of the whole galaxy 
population which met the above conditions, specifically 
elliptical or early type galaxies. Even in such objects, 
the presence of e.g. embedded disks caused the models 
to fail and produce residual maps marked by ubiquitous 
crosses or artificial features. The literature is literally 


full of such examples (e.g., Rothberg & Joseph 2004 

Reda et al.||20041 |Reda et al.||2005l IJanowiecki et al. 

20101 Peng et al.||20II| 

, IChies-Santos et al.||20II||Alamo- 

Martinez et al.|20I3[ G 

uerou et al.|20I5[ etc). Moreover, 


some of the correlations (or lack thereof) have been af¬ 
fected by the use of a formalism which fails to adequately 
capture the true isophotal shapes. 


In this work a new isophote fitting formalism is 
introduced, which is capable of modelling galaxies with 
significantly more complex or exotic isophote structures, 
with a particular focus on disk galaxies viewed edge-on. 
The paper is structured as follows. Section provides 
a short description of how deviations from perfect 
ellipticity are commonly expressed and modelled in 
the popular isophote analysis tool Ellipse. In Section 
a new angular metric is introduced, which better 
expresses isophote shapes in general and is particularly 
powerful in modelling disk galaxies viewed edge-on, as 
well as galaxies with X-shaped/peanut bulges. The 
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formalism is implemented and further demonstrated on 
a representative such galaxy. Section demonstrates 
the applicability of this formalism to any isophote 
shape, and explores scientific case-studies where this 
new technique can provide insight. In this section, four 
additional galaxies are modelled and the possibility to 
quantify the peanut/X-shape bulge of some galaxies 
from their photometric structure is briefly explored. 
Finally, in Section the paper reiterates its main 
conclusions and proposes further potential scientific 
applications of the method. 


2. ISOPHOTE FITTING 

Among the early efforts to express departures from 
p ure ellip t ical s hapes in mathematical form, the work 
of Carter (1978) was to become key in modelling more 
rehned and realistic galaxy isophotes. This work pro¬ 
posed adding perturbations to an ellipse as a function of 
azimuthal angle 0, in the same way as in a Fourier series 
decomposition 


= {ha) + X [Asin(n(/)) + S„cos(n(/))], (1) 


where 1(0) is the intensity profile along the isophote, 
expressed as a function of (central) azimuthal angle, 
{Iell) is the average intensity of the purely elliptical 
path, and the sum represents Fourier harmonic pertur¬ 
bations to (/ezz), with n being the harmonic (integer) 
order. Note that perturbing the intensity distribution in 
this way is equivalent to distorting the physical shape of 
the isophote. The two cases are used interchangeably 
throughout, for illustrative purposes. 

This formalism is particularly elegant as the coeffi¬ 
cients of the Fourier harmonics {A^ and BA carry phys¬ 
ical meaning. This is illustrated in Figure [l] for the first 
four harmonic orders. Here we observe how the harmon¬ 
ics bring corrections to the ellipse’s centre (n = 1), the 
ellipticity and position angle (n = 2), as well as captur¬ 
ing various types of asymmetries and the boxy or disky 
feature (the B 4 coefficient). 

These principles were later detailed into a full isophote 
fitting algorithn i in the seminal pap er of |Jedrzejewski| 


( 1987a[ ) (also in Jedrzejewski 


1987b). The reader is re- 


ferred to these works for a detailed description of the 
algorithm, which is only briefly outlined below. In this 
method, nested isophotes (not necessarily concentric) are 
fitted at several pre-defined points along the semi-major 
axis (denoted by a) of a galaxy in a CCD image. At 
each of these points, the isophote starts as a pure ellipse 
defined by initial guess-values for its geometric parame¬ 
ters, namely centre position {xQ^yo), position angle (PA) 
and ellipticity (e = 1 — 6/a, where b is the semi-minor 
axis). The image is then sampled along this elliptical 
path, giving the 1-dimensional intensity distribution as 
a function of azimuthal angle Idata{(t>)i which is first mod¬ 
elled by the right-hand side of Equation but restricted 
to n G {1, 2}, i.e.. 
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Figure 1. The physical significance of the first four harmonic 
corrections to an elliptical isophote. Positive coefficients are illus¬ 
trated in blue (thick) while negative coefficients in red (thin). The 
reference isophote is plotted in black (dotted) in each panel. 


After a least-squares minimisation of the quantity S 
in Equation the best-fit values of the harmonic coef¬ 
ficients are used to update the geometric parameters of 
the ellipse, and then {leii)- 

5 = X \hata{h) - I{h)f ■ (3) 

i 

The minimisation iterates until a minimum root mean 
square criterion is met. By only considering the first two 
harmonics (centre, PA and e), this section of the algo¬ 
rithm essentially computes the best-fiting pure ellipse on 
the data isophote corresponding to this semi-major axis 
radius. Einally, higher order harmonic perturbations to 
the ellipse (n > 3) are fitted (sequentially) through the 
same iterative minimisation. 

A powerful and robust implementation of this algo¬ 
rithm is in the IRAE (Image Reduction and Analysis Ea- 
cilit30) external software package Isophote, developed 
by the Space Telescope Science Data Analysis Systenj^ 
(STSDAS). The main tasks in this package (frequently 
referred to in the remainder of this work) are the isophote 
fitting task Ellipse, and the task Bmodel which builds 
a 2-dimensional galaxy model based on the parameters 


2 

h^) = (ha) + X Pnsin(n(/)) + BnCos{n(p)]. (2) 

n=l 


^ http://iraf.noao.edu/ 

^ http://www.stsci.edu/institute/software_hardware/stsdas 
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obtained by Ellipse. This software has had consid¬ 
erable success since its release, performing admirably 
when modelling the type of galaxies it was designed 
to work with, which are objcts characterised in general 
by relatively low ellipticity (e < 0.5) or, for higher e, 
with isophotes very close to elliptical (low harmonic am¬ 
plitudes). These are typically early-type galaxies with 
no (or with rather face-on) disks. The Ellipse task 
also performs quite well when modelling (late-type) disk 
galaxies viewed face-on, after the appropriate masking of 
biasing features such as dust lanes. 

However, for galaxies having more complex shapes or, 
in particular, for disk galaxies viewed at relatively high 
inclination or edge-on. Ellipse (and Bmodel) begin to 
break down and give rise to residual images marked by 
ubiquitous cross-like patterns or characteristic, alter nat- 


ing regions of excess and deficit light (see Sections 3.3 
and 

Obviously it is desirable to correctly model a galaxy’s 
light distribution. A good model provides physical 
insight into the galaxy itself, but moreover, subtracting 
an accurate model from an image gives the possibility 
to perform meaningful studies of substructures still 
remaining after the subtraction, such as star clusters, 
globular clusters or optical counterparts to X-ray/radio 
sources. In the case of edge-on galaxies, a model where 
the light intensity is significantly underestimated along 
the major axis of the galaxy introduces errors in the sur¬ 
face brightness profile, which is typically measured (and 
in fact output by Ellipse) along a, and is one of the 
main quantities of interest provided by isophote fitting 
programs. The surface brightness profile is essential in 
galaxy decomposition and is routinely used to quantify 
galaxy comp onents and ch aracterise structural scaling 
relations (see Graham||2013 review article and references 
therein). 



Figure 2. Panel a): Equal polar opening angles ( 0 ) define equal 
arc lengths on a circle (ci =02) but not on an ellipse (ei / 62). 
Panel b): When the same opening angles are transformed to the 
‘eccentric anomaly’ {'ip), they change as a function of (p such that 
they define equal arc lengths on an ellipse = ^2)^ longer 

on the circle {c'^ 7^ ^2). b therefore samples an ellipse uniformly. 


anomaly is frequently used in celestial mechanics to de¬ 
scribe the position of a point orbiting on an elliptical 
path, in canonical form (i.e where the angle is defined 
from the orbit centre, not from one focus), or more gen¬ 
erally to express the parametric equation of an ellipse. It 
is the ‘natural’ angular co-ordinate for ellipses, and the 7 /^ 
co-ordinate of a point on an ellipse is in fact related to its 
azimuthal/plane-polar angular co-ordinate (denoted as (j) 
throughout this paper) through the ellipticity e = 1—6/a, 
such that 


V, = -arctan ("5^) . (4) 

By transforming from 0 to t/, the isophotes are sam¬ 
pled uniformly along the entire azimuthal range (Figure 
right-hand panel), and the Fourier corrections are ex¬ 
pressed as 


3. A NATURAL ANGULAR METRIG 

The reason why isophote-fitting algorithms like 
Ellipse fail in the regime of high e and high harmonic 
amplitudes is traced back to the incorrect manner in 
which a quasi-elliptical isophote path is parameterised 
and sampled. The task Ellipse splits the isophote into 
sectors as it would a circle: it dividies the azimuthal 
range {(j) G [0;27r]) uniformly, in equal bins. While the 
sectors (arc lengths) corresponding to equal steps in (j) 
are all equal for a circle, they are not for an ellipse, but 
are actually longest along the majo^r axis and shortest 
along the minor axis (see Figure® left panel). This 
causes the isophote to be too coarsely sampled along 
a (in the case of edge-on galaxies, a corresponds to 
the plane of the disk, and thus needs to be very finely 
sampled as it contains most of the light), which then 
leads to bright streaks in residual images. 

3.1. The Eccentric Anomaly 

This work introduces a new isophote fitting formal¬ 
ism implemented in a new IRAF task ISOFIT, where the 
Fourier harmonics that quantify deviations from perfect 
ellipses are expressed as a function of an angular param¬ 
eter more natural to elliptical shapes, namely the ‘eccen¬ 
tric anomaly’, henceforth denoted by ip. The eccentric 


-f(^) = (leii) + X [Asin(nt/)) + B„cos(n?/))]. (5) 

n 

Furthermore, the Fourier harmonics correct the isopho- 
tal shape quite differently. This is illustrated in Fig¬ 
ure where the top panel shows how the cosine form 
of 7 / (which enters in Equation changes as the ellip- 
ticity changes from e = 0 ( 7 / = 0) to high ellipticity. 
The bottom panel illustrates how the cosine part of the 
n = 4 harmonic (i.e., boxyness/diskyness) acts sharply 
at the major axis (defined to correspond to 7 / = 0, 27r and 
ip = ir) and more softly along the minor axis {ip = 7r/2 
and pj = 37r/2), as the ellipticity increases. The more 
flattened (i.e., high e) an isophote, the sharper the har¬ 
monic correction around a and the softer the correction 
around 6. This is precisely what is required, since a disk 
galaxy viewed edge-on has a very bright but thin disk 
dominating the light along a, while along and around 
the minor axis 6 the light is dominated by a spheroidal 
bulge. The lower panel in Figure thus illustrates how 
the ‘lemon-shaped’ contours in the lower-right panel of 
Figurecan be made more ‘pointy’ along the major axis. 

A concrete example of harmonic terms characterising 
an isophote shape is given in Figure which also 
highlights the difference between using the generic polar 
angle (Ellipse) and eccentric anomaly (Isofit) to 
express the Fourier corrections. Here we notice how the 
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(p [rad./ 7i] 

Figure 3. Top: the cosine of the ‘eccentric anomaly’ '0, for differ¬ 
ent values of ellipticity (e) plotted against the polar angle 0. The 
amplitude of this harmonic {Bi in Figure Q controls the centre 
position of the isophote. Bottom: The cosine of 4'0, for the same 
values of e as in the top panel. The amplitude of this harmonic 
(B 4 in figure Q controls the boxyness/diskyness of the isophote. 

angle transformation (0 ^ tp) makes the A 2 correction 
correspond more to a rotation (change in position angle), 
and the B 4 correction resemble more a disky, lemon- or 
diamond-shaped contour. The lower two panels show 
a higher (n = 6) harmonic, which may carry physical 
meaning in edge-on galaxies with so-called X-shaped or 
‘peanut’ b ulge s. This tentative point is further explored 
in Section [4T| by modelling such a galaxy. 


3.2. Implementation in IRAF 

In order to model galaxy isophotes according to the for¬ 
malism discussed in the previous section, the source code 
of the IRAF package Isophote was modified at the har¬ 
monic fitting level (Ellipse) and 2D galaxy model build¬ 
ing (Bmodel). An angular parameter similar to the ec¬ 
centric anomaly has been used i n other galaxy-mo delling 
codes, most notably in Galfit (Peng et al. [2Q1Q| ), where 
it takes a slightly different forni (the Fourier formalism 
distorts the co-ordinate grid rather than the intensity 
distribution) and the authors consider only the cosine of 
the angle plus a phase in expressing the harmonics. 

The Ellipse task was modified and renamed Isofit in 
order to differentiate the new formalism from the stan¬ 
dard Ellipse. When fitting for the harmonic coeffi¬ 
cients, Isofit now defines the isophotes as a function of 
the eccentric anomaly tp by transforming the azimuthal 
angle (p according to Equation and applying the ap¬ 
propriate sign correction to each angle quadrant. Ad¬ 
ditionally, all harmonic correction^ of n > 2 are now 
added simultaneously^ in contrast to Ellipse, which fit¬ 
ted n > 2 harmonics sequentially. By fitting the har¬ 
monics simultaneously, all orders contribute together to 


^ In practce, n ~ 8 — 12 is the useful limit 



Figure 4. Examples of three harmonic corrections to an elliptical 
isophote, expressed as a function of polar angle (f) (left panels) and 
in the new formalism (as a function of ' 0 ; right panels), for values 
of e = 0.65 and A 2 = B 4 = Bq = 0.15, chosen to be relatively high 
to clearly illustrate the effect. In every panel, the reference ellipse 
is displayed with a dashed line. 


construct the isophote (just as in a Fourier decomposi¬ 
tion), rather than each individual order being a ‘best- 
attempt’ at capturing the dominant components of the 
galaxy light distribution. For example, in the particular 
case when the dominant component is a bright, thin disk 
along the major axis (i.e., an ecfee-on galaxy), we have 
seen in Figure and in Figure 0 (bottom two panels), 
and one can also deduce from Equation that every 
even-n harmonic provides a correction at the posi¬ 
tions corresponding to a {tp = 0 and tp = ir). Therefore, 
by fitting them sequentially each coefficient is optimised 
to best capture the disk feature rather than all orders 
contributing together to model less obvious features over 
the whole azimuthal range. 

In Ellipse, the harmonics are chosen by the user and 
their number is technically unlimited. This same flexibil¬ 
ity in terms of choice is kept in Isofit, though the simul¬ 
taneous fitting does limit the number of useful harmon¬ 
ics. Nevertheless, in practice the harmonic amplitude de¬ 
creases with increasing n and usually asymptotes to zero 
for n ^ 10 — 12 and above, even for extremely thin edge- 
on disks. Furthermore, in the majority of cases, odd-n 
harmonics bring little-to-no refinement to the isophote 
shapes (unless the galaxy has significant asymmetries), 
so by restricting the choice to only even orders, one can 
reliably fit harmonics up to n = 12 a nd ob tain accurate 
models, as is demonstrated in Section below. 

In addition to Ellipse, the auxilliary task Bmodel 
(‘build model’) was modified and renamed Cmodel 
(‘construct model’). The task Cmodel considers all har¬ 
monics included by the user (B MODEL only used n < 4) 
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and generates a 2D model of the galaxy based on the fit¬ 
ting results output by ISOFIT. The harmonic corrections 
are applied as a function of which is obtained as be¬ 
fore. In the regime e ^ 0, the standard Ellipse result 
is recovered, since in this limit pj = p. 


3.3. Detailed Example in Praetiee: the Edge-on Galaxy 
ESO 243-49 

The ideal test-cases, which best highlight the improve¬ 
ments brought by the formalism discussed in this work, 
are edge-on disk galaxies. These galaxies are charac¬ 
terised by diamond-like or lemon-like isophote shapes 
which deviate strongly from pure ellipses. They are the 
objects for which the standard version of Ellipse gives 
the least accurate results in both ID (systematically un¬ 
derestimated major-axis surface brightness profiles) and 
2D (cross-like features in residual images). 

The SO galaxy ESO 243-49 was chosen as the test-case 
for the new tasks IsOFiT and Cmodel. This edge-on 
galaxy has been intensly studied following the discovery 
of a source of hyperluminous X-ray emission (referred to 
as HLX-1) associated with a potential intermediate-mass 
black hole within the galaxy but outside the nuclear re¬ 
gion dFarrell et al.|2009 JWebb et al.|2010|). The fac t that 
HLX-1 was found to have an optical counterpart (Soria 
et al.||2QlQ ) makes this galaxy particularly interesting to 
model, as one of the aims of this work is to achieve high- 
quality 2D galaxy models and, after subtracting them 
from the image data, to perform accurate photometry 
on the remaining substructure (such as HLX-1) left be¬ 
hind in the residual image. Additionally, ES0243-49 also 
displays a relatively large number of star clusters both in 
the plane of the disk and outside, as well as a dust disk 
in the nucleus, also viewed edge-on. The image data 
for ESO 243-49 was retrieved from the Hubble Legacy 
Archiv^ and consisted of an iL-band image taken with 
the Wide Field Camera 3 (WFC3, IR channel, F160W 
filter, PI: S. A. Farrell) on the Hubble Spaee Teleseope 
(HST). 

To fully showcase the improvements brought by 
the new formalism, the image was modelled firstly 
with the standard Ellipse/Bmodel and then with 
Isofit/Cmodel. Both versions were run, where appli¬ 
cable, with identical initial conditions. 


3.3.1. Constructing the Isophotes 

An important observation concerning the construction 
of each isophote is that, due to the use of the polar an¬ 
gle (j) instead of pj to define the isophote in Ellipse, 
the high-order (n > 4) harmonics do very little to im¬ 
prove its shape (see Figured bottom-left panel). This is 
in part why, generally, higher order harmonics have not 
been used in the literature and it has been assumed that 
they are not important. This is illustrated in the top 
panel of Figure which shows a single isophote and its 
constituent even Fourier contributions for 2 > n > 10. 
This isophote was constructed from the harmonic coef¬ 
ficients computed by Ellipse while modelling ESO 243- 
49 The bottom panel in Figure displays an isophote 
constructed at the same semi-major axis but with the 

http://hla.stsci.edu 
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Figure 5. An example of how an isophote of ESO 243-49 (thick 
black line) is constructed from several harmonic orders, with the 
standard formalism (top) and the new /(b) formalism (bot¬ 
tom). The different harmonic orders, denoted by n, refer to Equa- 
tion|^(top) and Equation]^ (bottom). 

new formalism of IsOFiT proposed in this work. Here we 
see how each upper harmonic order brings a new level 
of refinement to the final shape, as one expects in fact 
from a Fourier analysis. The ability of each of these cor¬ 
rections to describe the actual isophote shape (at this 
semi-m ajor ax is) in ESO 243-49 will become apparent in 
Section 13.3.31 


3.3.2. ID Improvements 

A very useful quantity routinely extracted from galaxy 
images with the use of isophote fitting tools is, of course, 
the surface brightness profile /i(R) - the radial distribu¬ 
tion of light in the galaxy. This quantity is frequently 
used when performing one-dimmensional structural de¬ 
compositions (e.g., bulge/disk) of galaxies. Often times 
/i is measured along the major axis of the galaxy. For 
edge-on galaxies, this is precisely the locus where light 
is systematically underestimated when using improperly 
defined isophote shapes. In such cases (and the test- 
case ESO 243-49 is a prime example) one must resort to 
estimating ja{R) by taking a crude ‘cut’ along the im¬ 
age major axis, which is undesirable for several reasons. 
Firstly, this method is noisy, because of e.g., dust, star 
clusters or foreground stars which happen to be located 
in the plane of the disk. Further, when taking a cut, 
one obtains a single measurement of the brightness at 
each radial position, whereas by fitting the isophote, one 
obtains a value for the brightness at a given R which 
is the average value across the whole azimuthal range. 
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Secondly, a cut is a less accurate representation of the 
major axis surface brightness profile because, in general, 
the position angle of the isophotes does not remain con¬ 
stant with increasing radius: real isophotes often rotate 
with in crea sing i?, following features such as bars (see 
Section so a ID cut capturing this effect must be a 
curve, not a line. 

Figure shows how the use of the old formalism in 
Ellipse systematically underestimates the major axis 
surface brightness profile of the test-case galaxy ESO 
243-49, compared with a (constant PA) direct cut taken 
along the disk plane and assumed to be a more realistic 
approximation of the true /i(i?maj)- The discrepancy is 
most prominent across the range (5 < i^maj ^ 30 arc- 
sec) where the disk light dominates. The formalism pre¬ 
sented in this work does remarkably better and is the 
closest and most accurate representation of for 

the reasons explained above. The discrepancy between 
the profiles computed with Ellipse and Isofit is max¬ 
imal at i^maj = lb.5 arcsec, where Ellipse underesti¬ 
mates /i by 0.77 mag arcsec”^, which is a factor of two 
in surface brightness. 

Such a discrepancy propagates into the bulge/disk de¬ 
compositions and leads to wrong results. Judging by the 
shape of the underestimated profile (green dashed curve 
in Figure]^, the disk component appears to follow an ex¬ 
ponential form. If modelled as such, it would eat into and 
impact the fit of the central bulge (its brightness, concen¬ 
tration, effective radius etc.). When considering the cor¬ 
rect profile (black solid curve in Figure]^, the disk com¬ 
ponent is in reality more appropriately described by an 
inclined disk model (a Sersic function with nsersic < 1; 
Pastrav et al.||2013|), which leads to quite a different de¬ 
composition result, with a brighter and more concen¬ 
trated (higher Sersic index) bulge. When performing 
bulge/disk decompositions, it is common for the major- 
axis surface brightness profile, /i(i?maj): to be mapped 
onto the ‘equivalent’ axis i?eq- This mapping is equiv¬ 
alent to transforming each isophote of semi-major axis 
i?maj into a circle of equal enclosed area, which then has 
a radius i?eq- This way, one can use azimuthal symmetry 
to compute integrated quantities from /i(i?eq), including 
the magnitudes of the model components that are fit to 
this profile. This is fully consistent with the formalism 
implemented in this work as the Fourier harmonic per¬ 
turbations change the shape of the isophote but conserve 
its area. This is demonstrated in the Appendix. 

A final note on the ID aspect of this work is that 
using the correct angular metric to express the isophotes 
makes it possible to obtain the surface brightness profile 
not just of the major axis, but also the minor axis or 
in fact any direction, directly by using (leii) and the 
harmonic coefficients (all provided by Isofit) in Equa¬ 
tion for any azimuthal angle. This was not possible 
before, and we can readily see it just by inspecting the 
single isophote shown in the upper panel of Figure 
(the shape as a function of azimuthal angle is wrong). 



Figure 6. Major axis surface brightness profile of ESO 243-49. 
The standard Ellipse task systematically underestimates the pro¬ 
file in the range 5 — 30 arcsec (where the disk dominates), by as 
much as a factor of 2 in surface brightness occurring at Rmaj = 16.5 
arcsec. See the lower panels of Figure to understand this. 


proper angular metric employed, the residuals obtained 
from subtracting these 2D models from galaxy images 
are notorious for displaying obvious artificial features, in 
particular for systems with strong departures from pure 
elliptical shapes. In this subsection Bmodel and the new 
task Cmodel, are both run on the test-case galaxy ESO 
243-49. The input information for Bmodel consisted 
of the best-fit generated by Ellipse, while the input for 
Cmodel consisted of the best-fit generated by Isofit. 
The results are shown in Figure 

We can clearly recognise in the upper-left panel of Fig¬ 
ure the wrong shape of a disky isophote resulting from 
imprOTerly defined ellipses (we noticed this in Figures]^ 
andl^. The inevitable consequence is that when such a 
mod3 is subtracted from the galaxy image, the resulting 
residual map is contaminated by artificial features (for 
ESO 243-49, the residual image is dominated by such 
an artificial feature along the disk plane, which is quite 
obvious in the bottom-left panel of Figure [^. 

It is immediately apparent that Cmodel (top-right 
and bottom-right panels in Figure achieves a re¬ 
markably better representation of me galaxy light 
distribution than the previous standard (Bmodel). 
With Cmodel it is now possible to construct realistic 
and accurate 2D models of galaxies. Having such models 
can bring considerable insight into the structure and 
morphology of galaxies and also enables meaningful 
studies of the substructure left behind after subtract¬ 
ing the accurate model from the image. As to the 
latter, Cmodel makes it possible to easily identify 
the optical counterpart of the X-ray source HLX-1 and 
compute its near-infrared brightness. This wa s done 
with the so ftware Aperture Photometry Tool ( |Laher| 


et al. 2012), resulting in the apparent iJ-band magni- 
tude tuh = 22.68±0.27 (in the Vega magnitude system). 


3.3.3. 2D Improvements 

The auxilliary task Bmodel in the IRAF package 
Isophote has been routinely used in the past to gen¬ 
erate two-dimensional models of galaxies, based on the 
ID information computed by Ellipse. Due to the im- 


4. NEW SCIENCE - CASE STUDIES 

This section demonstrates the much broader applica¬ 
bility of the new technique to describing virtually any 
galaxy. Four additional galaxies are modelled, chosen to 
be representative of a specific class that presents scien- 
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Figure 7. Top 3 panels: 2D model of the galaxy ESO 243-49 constructed with old formalism of Ellipse and Bmodel (top-left); HST 
image of ESO 243-49 (top-centre); 2D model constructed with ISOFIT and Cmodel (top-right). Bottom 2 panels: Residual images obtained 
after subtracting the BMODEL-based galaxy model (bottom-left) and the new CMODEL-based galaxy model (bottom-right) from the galaxy 
image. The position of the optical counterpart of HLX-1 is marked with the cyan circle in both residual maps. The higher harmonic terms 
were turned on and used to construct both models. However, the old formalism used in Ellipse/Bmodel fails due to the reasons detailed 
in Section [ 3 ] 


tific interest and is the object of active current research. 
The following is by no means an exhaustive study of 
each class, but it is meant to merely illustrate the use¬ 
fulness of the method introduced in this work, and to 
highlight some tentative ways to investigate such galax¬ 
ies and quantify their properties. 


4.1. Peanut/X-shaped Bulges & Boxy+Disky Isophotes 

The bulges of many nearby galaxies have been ob¬ 
served to display non-elliptical shapes, particularly ob¬ 
vious in (but not restricted to) galaxies where the bulge 
is embedded in a large scale disk. These bulges are com¬ 
monly referred to as boxy, peanut-shaped or X-shaped, 
and have been amply studied both observa t ionally and 
with simulations (e.g de Vauconlenrs| 1974[ Combes fc 
Sandersl [19811 IBnreau hreeman||1999[ [Bureau et al. 


2006 


- see 


Laurikainen &: Salo||2015| for a recent review). 


In terms of formation scenarios, it is believed that this 
structural feature has a kinematical origin, resulting from 
the re-organisation of stellar orbits within the bulge, and 
is possibly related to the formation of bars. Because 
peanut/X-shaped bulges are especially obvious when the 


galaxies are viewed edge-on, these objects make ideal 
candidates for the modelling technique described in this 
work. 

A representative galaxy, known to have an X-shaped 
bulge, was chosen to be modellled in this section, namely 
NGC 2549. The imaging data was retrieved from the 
HLA, and consisted of an i?-band HST image taken with 
the Wide Field and Planetary Camera 2 (WFPC 2 , filter 
F702W; PI: W. Jaffe). As before, the galaxy was mod¬ 
elled with the old method (standard Ellipse , Bmodel) 
and the new technique (ISOFiT, Cmodel), both run with 
identical starting conditions. The results are displayed 
in Figure 

As mentioned in Section [^ one of the goals of this work 
is to investigate whether a galaxy’s more complex mor¬ 
phological features are encoded in its isophotal structure, 
and can therefore be quantified through the Fourier co¬ 
efficients. While the H 4 harmonic coefficient has been 
shown in the past to correlate with structural and phys¬ 
ical properties of elliptical galaxies, the higher orders 
and late-type galaxies have been left quite unexplored 
(largely due to the improper expression of the isophote 
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Figure 8. NGC 2549: a galaxy with a known X-shaped bulge. The panels are arranged analogously to FigureThe upper panels use a 
segmented colormap to better highlight the X-shaped bulge feature, and only 3 contours are overplotted: the inner and outer black contours 
have boxy and disky shapes, respectively. The middle contour (green) is a combination of both and needs higher orders to be properly 
modelled. 


contours). With this new formalism it is now possible to 
explore this parameter space. 

If we observe the isophote contours of NGC 2549 in the 
upper panels (centre and right) of Figure we notice 
how the isophotes transition from boxy (inner) to disky 
(outer). The shapes at the transition are not elliptical 
however (even though ^4 = 0 ), but rather a combina¬ 
tion of boxy and disky (green contour), reminiscent of 
the ^e-corrected ellipse illustrated in Figure (bottom 
right panel). The transition occurs precisely in the radial 
range where the X-shape/peanut feature is most obvious 
in the image. Since it is clear that higher-order harmon¬ 
ics are necessary to capture this morphological feature, 
the harmonic coefficients actually ‘measure’ the peanut- 
ness in this way. This is an important result: with the 
degree of peanutness quantifiable this way, it is now pos¬ 
sible to directly and quantitatively compare ob servations 


of real galaxies with simulation s (such as e.g., Athanas- 
sonla|[2Q^ , Athanassonla||2Ql^ . 


A second galaxy, NGC 361U, which displays a combi¬ 
nation of boxy and disky isophotes was modelled in this 
section. This galaxy is consistent with having a thick 
boxy bulge, i.e., the extent of the disk component is less 
than or comparable to that of the bulge. In this work 


the terminology ‘boxy/peanut (B/P) bulge’ is avoided, 
however, precisely because boxyness is expressed through 
^ 4 , whereas the isophotes of these structures are associ¬ 
ated with different/additional harmonics. In stead, these 
structures are referred to as X-shaped/peanut (X/P) 
bulges. The data was again retrieved from the HLA, 
and consisted of an /-band image taken with WFPC2 
(filter F814W; PI: B. Whitmore). The resultind models 
and residual maps corresponding to this galaxy are dis¬ 
played in Figure Once again, standard Ellipse can 
not model the higher harmonic orders, and the residual 
image is dominated by an artificial streak along the disk 
plane. 

Figure 


10 displays the B 4 and Bq coefficients of NGG 


2549 and NGG 3610, computed with ISOFIT, as they 
change with radius (along the major axis). The coef¬ 
ficients are re-normalised by the local gradient and local 
semi-major axis length (Equation |^, such that 


Bn —Bn ( ci-p;— 
oa 


-1 


(6) 


where a is the local semi-major axis length of the 
ellipse. This is a common way of expressing the coeffi- 
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Figure 9. NGC 3610. The panels are arranged analogously to Figure [7| The isophote shapes require higher-order harmonics to properly 
capture the combination of both boxyness (‘square’ along the minor axisjand diskyness (‘pointy’ along the major axis due to the embedded, 
nearly edge-on disk) within the same isophote. 




Figure 10. The fourth and sixth Fourier coefficient profiles corresponding to the is^hotes of NGC 2549 (left panel) and NGC 3610 (right 
panel), normalised by the local semi-major axis and intensity gradient (see Equationand plotted as a function of semi-major axis length. 
In NGC 2549 there is a transition between boxy (negative B 4 at Rmaj $ 50 pixels) and disky (positive B 4 at Rmaj ^ 50 pixels) isophotal 
shapes, and the higher order Bq term is required to capture the ‘peanut’ feature. NGC 3610 has a similar transition in B4 and also requires 
a Bq contribution. 
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cients as (dimensionless) deviations from ellipticity, as 
opposed to concrete intensity corrections. The transition 
between boxy and disky isophotes {B 4 crossing 0) occurs 
in both galaxies, and also Bq has non-zero amplitude for 
most of the radial range, which indicates that it plays 
a part in capturing the isophote shape. Such accurate 
information was not available with the old formalism 
used in IRAF. 


4.2. Photometry of Galaxy Substructure 

Perhaps the most basic use of a high-quality model of a 
galaxy is the study of residual substructure, such as star 
clusters, globular clusters (GCs), satellites, streams etc. 
Once the model is subtracted from the original image, 
all of this substructure is lef t beh ind in the residual im¬ 
age. We have seen in Section [3~3| that Cmodel makes it 
possible to perform photometry on residual images quite 
straightforwardly with tools such as the Aperture Pho¬ 
tometry Tool. In this section we qualitatively explore 
this further, and extend the study to modelling galax¬ 
ies with more complex isophotal structures such as mis¬ 
aligned bars. To begin with, the interesting barred lentic¬ 
ular (SBO) galaxy NGC 936 was the first candidate to be 
modelled. Due to this galaxy’s orientation (slightly in¬ 
clined but close to f ace-ou view) and compo nents (bulge. 


disk and barlens; see 

Laurikainen et al.| 

2011), its isophote 

structure is relatively complicated. ' 

I'he data was re- 


trieved from the publicly available archive of the Spitzer 
Survey of Stellar Structure in Galaxies (S4GJ^, and con¬ 
sisted of an infrared (3.6/i) Spitzer observation (IRAC 
1 instrument, PI: K. Sheth). Figure pT] displays the re¬ 
sults: the isophote shape starts off close to elliptical in 
the central regions, then it transitions to lemon-shaped 
and parallelogram-shaped (when capturing the barlens) 
and finally returns to a shape close to elliptical, when 
the barlens no longer dominates and the disk takes over. 
The isophotes also rotate in position angle with increas¬ 
ing radius, a process which is again driven by the bar¬ 
lens. Because of all these strong deviations from elliptic- 
ity, the 2D model built with Bmodel fails to reproduce 
the light distribution in this galaxy and after its sub¬ 
traction from the image, a residual map remains that 
is heavily contaminated by twisting features, both pos¬ 
itive and negative. The much more realistic 2D model 
built with Cmodel (top right) and the associate resid¬ 
ual map (bottom right) demonstrate that the formalism 
proposed here is not limited to edge-on systems, but can 
competently handle any type of galaxy. With such good 
quality results, it is now possible to make stronger claims 
about galaxy morphology and also gain better insight 
into the substructure than before. The previous contam¬ 
inated residual images made it difficult to analyse globu¬ 
lar clusters because the artificial features often hid them 
in oversubtracted regions or covered them through the 
bright artificial streaks. This not being a problem with 
our new approach, it is now possible to perform GC num¬ 
ber statistics (e.g., GG number density radial profiles) as 
well as more accurate photometry. 

The second galaxy modelled in this section was 
deliberately chosen from a survey dedicated to GC 


studi e s, namely the SLU GG^ survey (Usher et al. 


2012, [Brodie et al. 2012). Tne galaxy m question 


NGC 4111, is another example of a disk viewed at close 
to 90° inclination, and is in fact the ‘flattest’ galaxy 
modelled in this paper. The data was again retrieved 
from S4G archive, and consisted of an infrared (4.5/i) 
Spitzer observation (IRAG 2 instrument, PI: K. Sheth). 
The results are displayed in Figure where again we 
notice the contaminated residual image produced with 
Bmodel as well as the ‘cleaner’ residual obtained with 
with Cmodel. 


5. SUMMARY AND OUTLOOK 

This paper proposes a new formalism in which the 
isophotes of galaxies are expressed as a function of an 
angular parameter more natural to ellipses, namely the 
‘eccentric anomaly’. This replaces use of the polar an¬ 
gle and significantly improves on past algorithms, such 
as the IRAF tasks Ellipse and Bmodel. Isohotes with 
significant deviations from pure ellipses can now be ac¬ 
curately modelled, opening the door to new science. 

This new formalism is implemented in the tasks 
Ellipse and Bmodel, which are renamed Isoeit and 
Cmodel, respectively, in order to differentiate the new 
implementation from the old standard. The method is 
tested on an edge-on galaxy, ESO 243-49, and is found to 
bring considerable improvements in both ID (major axis 
surface brightness profile) and 2D (model and residual 
image). Specifically, with ISOEIT the true surface bright¬ 
ness profile of ESO 243-49 is recovered, in shape and 
normalisation, while with Ellipse it was systematically 
underestimated, the discrepancy reaching a maximum 
of 0.77 mag arcsec“^ (a factor of ^ 2 in brightness) at 
Rmaj = 16.5 arcsec for the galaxy in question. In addi¬ 
tion, the second new task Cmodel was found to produce 
a realistic 2D representation of the light distribution, 
which, after being subtracted from the image, left behind 
a high quality residual image on which photometry was 
possible. The iL-band magnitude of the optical counter¬ 
part of HLX-1 was measured to be mn = 22.68 ± 0.27 
mag. 

This work makes it possible to study a wide range 
of photometrically interesting objects, and two partic¬ 
ular scienctific applications (beyond galaxies with near 
edge-on disks) were further considered: galaxies with X- 
shaped/peanut bulges, and photometry of galaxy sub¬ 
structure. By expressing isophotes with a more natu¬ 
ral metric, the higher order Fourier moments now carry 
meaning and can potentially quantify morphological fea¬ 
tures such as the peanut feature in galaxy bulges. This 
conclusion, however, is only based on the two galaxies 
modelled in the paper and an exhaustive and more quan¬ 
titative study on a larger sample is reserved for future 
work to confirm this hypothesis. While the full power of 
this new formalism is most obvious in edge-on galaxies, 
it is not limited to these systems, and in fact is fully 
applicable to any galaxy. The task Isoeit is capable 
of capturing complex isophote shapes with high devia¬ 
tions from ellipticity and, in the limit of low ellipticity or 
low deviations from ellipticity, the results of Ellipse are 
fully recovered. 


^ http://irsa.ipac.caltech.edu/data/SPITZER/S4G/ 


® http://sluggs.swin.edu.au/ 
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Figure 11. NGC 936: an SBO galaxy with a barlens and complex isophotes. The panels are arranged analogously to Figure The 
barlens induces strong departures from pure elliptical shapes, as well as substatial (close to 90°) isophote rotation. 


The case-studies discussed in the above sections repre¬ 
sent a small sample of the much wider range of science 
made possible by the better mathematical description of 
isophotal structure used in ISOFIT. One such application 
involves the modelling of compact, low-mass early type 
galaxies. Current i nvestigations on such galaxies (e.g., 
Guerou et al. 2015) still model them with Ellipse and 
thus fail to capture interesting features such as embedded 
disks, which in turn leads to the familiar cross-pattern 
in residual maps. Because of this, often times the har¬ 
monics are deliberately not included, and the presence 
of disks/diskyness is pointed out by the cross pattern in 
the residual image. 

Additionally, recent studies have shown that the bulge 
of the Milky Way has an X/peanut-shaped structur e 


dWegg et al.|2Q15 , Di Matteo et al.|2Q"T5 Qin et al.|2Q15 ) 
Since this makes it the nearest X/F bulge in an edge-on 
galaxy, it is a prime candidate to study with IsOFiT and 
investigate potential connections between the X/P char¬ 
acteristic and n > 4 Fourier coefficients. 

Finally, this met hod can be applied to the interesting 
shell- galaxies (e.g., Forbes et al. 1994 Del Burgo et al. 
2008). These are elliptical galaxies marked with shell¬ 
like or ring-like structures, thought to be the remnants 
of past merger events. 


All of these constitute extensive projects in their own 
right and are deferred to future work. 
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Figure 12. NGC 4111: an edge-on galaxy surrounded by globular clusters. The panels are arranged analogously to Figure The 
artificial cross in the residual obtained from Bmodel (bottom left) makes it hard to detect all the globular clusters, and contaminates the 
computation of their brightness. Application of the new Cmodel obtains a much better residual image (bottom right). 
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APPENDIX 

The formalism proposed in this work is fully consistent with re-mapping the major-axis surface brightness profile 
/^(^maj), computed by ISOFIT, onto the ‘equivalent axis’ Req. The ‘equivalent radius’ i?eq is the radius of the 
‘circularised’ isophote of semi-major axis length i?maj such that circle encloses an equal area as that of the isophote. 
This transformation is routinely performed on /i(i?maj) following a bulge/disk decomposition of because it 

allows the use of circular symmetry in computing the integrated surface brightness (i.e., the magnitude) of the model 
components. The deviations from a purely elliptical shape brought by the IsOFiT formalism conserve the enclosed 
area, as is demonstrated below. 

I.) Consider a surface enclosed in a circle of radius a. In plane-polar co-ordinates, the surface element dS is given 
by: 


dS = rdrd(l). (A. I) 

The total area obtained by integrating dS in the range r G [0, a] and 0 G [0, 27r] is the usual 

pa p27t 

A= / rdrd0 = 7ra^. (A.2) 

Jo Jo 

II. ) Consider re-scaling the circle along the ^-axis by a factor h/a. The resulting shape is an ellipse with semi-major 
axis a and semi-minor axis b (Figure]^ panel a). Moreover, the area A is also re-scaled by the same factor: 
^eii = A bja = 710 “^ {bja) = 7ra6, the usual expression for the area of an ellipse. Further, a point F on the circle 
(defined py a central angle 0) becomes the point F' on the ellipse, and is now defined by a central angle ip. The two 
angles are related to each other through the axis ratio 6/a or ellipticity e = I — a/6, through Equation]^ This is in 
fact the geometric definition of the ‘eccentric anomaly’ 7 /. 

III. ) Consider again the circular surface, but perturbed with a set of sinusoidal waves, as a function of polar angle 
(e.g.. Figure]^ panel b), such that 


r(0) = r[I -h A^sin(n0) -1- 5nCos(n0)], (A-3) 

where are constant coefficients and n is an integer. The new surface element (denoted as dS') becomes: 

dS' = r[I -h A^sin(n0) + B^cos(n0)]drd0. (A-4) 

As before, the total area of the perturbed shape (denoted as A') is given by integrating the surface element in the 
same limits: 


n 27r 

r[I + A^sin(n0) + 5^cos(n0)] drd^ 

. 

pa p2'k pa p2'k 


pa pztt pa pztt pa p2'k 

/ / r drdcj) -h / / rsin(n0) drdcj) + / / rcos(n0) dvdcp 

Jo Jo Jo Jo Jo Jo 

p27t ^2 r27T 

/ sin(n(/))d(/) cos(n(/))d(/) 

Jo 2 Jq 


= na +- 


= TTO^ -h 7^[-I + 1+0-0] 

Zn 

.2 _ 


(A.5) 


= 7ra^ = M, V n 


The sinusoidal functions change the area as a function of p but the net change is zero over the full range [0,27r], i.e., 
the positive ‘bumps’ are exactly cancelled out by the negative ‘dips’ over the whole azimuthal range, independently 
of the choice of n. 


IV.) Consider now the perturbed shape ‘squeezed’ along the y direction by a factor 6/a (Figure]^ panels b and c), 
in the same way as the circular shape was re-scaled (panel a). A point P on the perturbed circle, wnich was described 
by the polar angle p (panel b), becomes the point P' on the perturbed ellipse (panel c), described by the eccentric 
anomaly as before. This is at the heart of the the new formalism - a Fourier wave on a circle needs a different 
co-ordinate system on the corresponding ellipse. As a consequence of the re-scaling (‘squeezing’), all areas (in panel 
b) are reduced by the same factor (6/a) in panel c. The area inside and outside of the ellipse is thus reduced by 6/a, 
but there is still, obviously, an equal amount of area inside and outside. The area of the ellipse is therefore conserved 
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Figure Al. Panel a): When the circle of radius a is ‘squeezed’ along the y—axis by a factor 6/a, it results into an ellipse of semi-minor 
axis 6 and semi-major axis a. The point F on the circle (defined by cf)) co rresp onds to F' on the ellipse, which is described by the ‘eccentric 
anomaly '0. Panel b): The circle in panel a) distorted through Equation |A.3[ with n = 4, An = 0, Bn =0.1 (thick line). A point P along 
the circumference is defined by the polar angle cf). Panel c): When the distorted circle in panel b) is ‘squeezed’ into a distorted ellipse, the 
same point (now denoted by P') is now defined by the ‘eccentric anomaly’ 'ijj. 


following the perturbation. 

V.) From the output of ISOFIT {{leii) , ^maj and An^Bn - see Sectionfor the explanation) the major axis surface 
brightness profile p{Rmaj) can be constructed in 2 equivalent ways: 

1. By applying the perturbations to the intensity/surface brightness (/^=o ^ {Rii)) ~ keeping the isophote shape 
elliptical and distorting the intensity distribution along its circumference. This is how IsOFiT works. 

2. By keeping (leii) fixed and applying the corrections to i^maj “ distorting the ellipse shape and keeping the 
intensity fixed. This step additionally requires the intensity gradient, which is also provided by IsOFiT (Ellipse 
also computes all these quantities, though with incorrect 5^), to re-normalise the coefficients (such that they 
correspond to i^maj perturbations). 

Performing the distortion the second way is analogous to case IV.): a given isophote is a distorted ellipse whose 
area is the same as that of the original, pure ellipse {Aiso = Aeii = ^rab = 7ra^(l — e) = — e)). The intensity 

associated with this isophote area is {leii) (which is converted into the surface brightness fieii)- Therefore, if the area 
is expressed as a circular area, its radius i?eq can be obtained as: 

= (A.6) 


-^eq -^maj ^ 


(A.7) 


With the surface brightness profile now expressed along the equivalent axis, it is straightforward to compute 
magnitudes in the usual fashion. 
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